arXiv:1505.04103vl [cs.NA] 15 May 2015 


A splitting scheme to solve an equation for fractional 
powers of elliptic operators^ 

Petr N. Vabishchevich^’^’* 

^Nuclear Safety Institute, Russian Academy of Sciences, 52, B. Tulskaya, 115191 Moscow, Russia 
^North-Eastern Federal University, 58, Belinskogo, 677000 Yakutsk, Russia 


Abstract 

An equation containing a fractional power of an elliptic operator of second or¬ 
der is studied for Dirichlet boundary conditions. Finite difference approximations in 
space are employed. The proposed numerical algorithm is based on solving an auxil¬ 
iary Cauchy problem for a pseudo-parabolic equation. Unconditionally stable vector 
additive schemes (splitting schemes) are constructed. Numerical results for a model 
problem in a rectangle calculated using the splitting with respect to spatial variables 
are presented. 
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1. Introduction 


Non-local applied mathematical models involving fractional derivatives in time and 
space are actively discussed HSIil at the present time. A lot of models in applied 
physics, biology, hydrology and finance include both sub-diffusion (fractional in time) 
and supper-diffusion (fractional in space) operators. Supper-diffusion problems are 
treated as evolutionary problems with a fractional power of an elliptic operator. 

For solving problems with fractional powers of elliptic operators, we can apply 
finite volume or finite element methods oriented to using arbitrary domains and irreg¬ 
ular computational grids lll2l|l4|] . A computational realization is associated with the 
implementation of the matrix function-vector multiplication. For such problems, vari¬ 
ous approaches iQl were developed. The application of Krylov subspace methods with 
the Lanczos approximation for solving systems of linear equations associated with the 
fractional elliptic equations is discussed in A comparative analysis of the con¬ 
tour integral method, the extended Krylov subspace method, and the preassigned poles 
and interpolation nodes method for solving space-fractional reaction-diffusion equa¬ 
tions is presented in fl. The simplest variant is based on the explicit construction of 
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the solution using the known eigenvalues and eigenfunctions of the elliptic operator 
with diagonalization of the corresponding matrix Sii. Unfortunately, all these 
approaches demonstrates high computational costs for multidimensional problems. 

We have proposed lE^I a numerical algorithm to solve an equation for fractional 
powers of elliptic operators that is based on a transition to a pseudo-parabolic equation. 
For an auxiliary Cauchy problem, the standard two-level schemes are applied. The 
computational algorithm is simple for practical use, robust, and applicable to solving 
a wide class of problems. A small number of time steps is required to find a solution. 
This computational algorithm for solving equations with fractional powers of operators 
is promising for transient problems. 

In the study of difference schemes for time-dependent problems of mathemati¬ 
cal physics, the general theory of stability (well-posedness) for operator-difference 
schemes |II1[I3] is in common use. At the present time, the exact (matching nec¬ 
essary and sufficient) conditions for stability are obtained for a wide class of two- and 
three-level difference schemes considered in finite-dimensional Hilbert spaces. We em¬ 
phasize a constructive nature of the general theory of stability for operator-difference 
schemes, where stability criteria are formulated in the form of operator inequalities, 
which are easy to verify. 

In numerical solving initial-boundary value problems for multidimensional PDFs, 
great attention is paid to the construction of additive schemes || 1 3ll2ill . The transition 
to a chain of simpler problems allows us to construct efficient difference schemes. Here 
we speak of splitting with respect to the spatial variables. In some cases, it is useful to 
separate subproblems of distinct nature. In this case, we say about splitting into phys¬ 
ical processes. Such schemes appear in the solution of unsteady problems for systems 
of interconnected equations. There are actively discussed regionally additive schemes 
(domain decomposition methods), which are focused on designing computational algo¬ 
rithms for parallel computers. Iterative methods for solving steady-state problems are 
often treated as pseudo-time evolution methods for solving time-dependent problems. 
Many iterative methods may be associated with the use of certain additive schemes. 

In this paper, we consider issues of constructing unconditionally stable schemes 
for approximate solving problems with fractional powers of elliptic operators on the 
basis of a pseudo-time evolutionary problem. The success is achieved through the 
use of vector additive difference schemes of multicomponent splitting US. The 
original problem is reformulated as a vector problem. In this case, instead of a single 
approximate solution, we search a vector of approximate solutions. The corresponding 
additive schemes are schemes of full approximation, i.e., at each time level, we search 
the approximate solution of the problem. 

The paper is organized as follows. The formulation of a steady-state problem for 
a space-fractional elliptic equation is given in Section|2] Finite difference approxima¬ 
tions in space and the standard two-level schemes to solve an auxiliary Cauchy problem 
for a pseudo-parabolic equation are discussed in Section[3 In SectionlH we construct 
a special additive difference scheme for time-stepping and investigate its stability. The 
results of numerical experiments are described in Section|5] 
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2. Problem formulation 

In a bounded polygonal domain Q c with the Lipschitz continuous boundary 
dn, we search the solution for a problem with a fractional power of an elliptic operator. 
Introduce the elliptic operator as 


= -div^(jc)grad M + c(x)u (1) 

with coefficients 0 < < k(x) < k 2 , c{x) > 0. The operator ^ is defined on the set of 

functions u(x) that satisfy on the boundary dQ. the following conditions: 

u{x) = 0, X e dQ.. (2) 


In the Hilbert space "K = l-qiQ), we define the scalar product and norm in the 
standard way: 

(u,w) - I u(x)w(x)dx, \\u\\ - (u,uy^^. 

Jn 

In the spectral problem 

^(fic = Aktpk, X eQ, 
tpk-O, a: 6 dQ, 

we have 

Ai < A 2 < ..., 

and the eigenfunctions ipk, WfkW - k- 1,2,... form a basis in L 2 (Q). Therefore, 


k=l 

Let the operator be defined in the following domain: 

CO 

D{.9l) = (m I u{x) e L2{Q), \(u, ipk)?'Ak < <»}. 


k=Q 


Under these conditions : L 2 {Q) —> L 2 (Q) and the operator Jl is self-adjoint and 
positive definite: 

> All, Ai > 0, (3) 

where I is the identity operator in 'K. In applications, the value of Ai is unknown (the 
spectral problem must be solved). Therefore, we suppose that 6 < Ai in (l3]l. Let us 
assume for the fractional power of the operator Jl: 

00 

= Yj('^^‘Pk)AkiPk. 

k=0 


More general and mathematically complete definition of fractional powers of elliptic 
operators is given in 11241] . The solution u{x) satisfies the equation 




(4) 


under the restriction 0 < cr < 1. 
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3. Numerical algorithm 


We consider the simplest case, where the computational domain Q is a rectangle: 


Q = {jc IJC = (xi,X2), 0 < Xk <lk, k - 1,2}. 

To solve approximately the problem (|4|i, we introduce in the domain Q a uniform grid 

(n={x\x = (xi,X2), Xk^khk, ik^0,l,-,Nk, Nkhk = h, k ^ 1,2], 

where aJ = cu U dot and a> is the set of interior nodes, whereas dot is the set of boundary 
nodes of the grid. For grid functions y(jc) such that y(jc) - 0, x i a>, we define the 
Hilbert space H - Li (u), where the scalar product and the norm are specified as 
follows: 

= ^y(x)w{x)hih2, llyll = ■ 

X&CJ 

For the discrete operator A, we use the additive representation 

2 

A-'^Ak, X eoj, 

k=\ 

where Ak, k - 1,2 are associated with the corresponding differential operator with the 
second derivative in one direction. 

For all grid points except adjacent to the boundary, the grid operator Ai can be 
written as 

Aiy = - -j-rkixi +0.5/ii,/i2)(y(xi + huh 2 ) -y(x)) 
h] 

+ -^kixx - Q.5hx,h2)(y{x) - yixx - hi,h 2 )) + c(x)y{x), 
h\ 

xeuj, x\+0.5hi, xi /j - 0.5/ji. 

At the points that are adjacent to the boundary, the approximation is constructed taking 
into account the boundary condition (|2|l: 

Aiy = - ^k{xi + 0.5hi,h2){y{xi + hi,h 2 ) - y{x)) 
h\ 

+ -^k{x\ - 0.5h\,h2)y{x) + c{x)y{x), 
h\ 

X e 10, x\ - 0.5hi, 


Aiy =-^k(xi + 0.5huh2)y(x) 

h\ 

1 

+ -^k(xi - 0.5huh2){y{x) - y(xi - /ri, /ra)) + c{x)y{x), 

h\ 

X € CO, xi — li - 0.5hi. 
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Similarly, we construct the grid operator A 2 . For the above grid operators (see, e.g., 
I Tsl 1 ^ ), we have 

Ak^Al>6kI, 6k >0, k^l,2. 


where I is the grid identity operator. Because of this, the discrete operator A is self- 
adjoint and positive dehnite in H: 


2 

A^A’-^dl, (5) 

k=l 

For problems with sufficiently smooth coefficients and tfie right-hand side, it approxi¬ 
mates the differential operator with the truncation error O (ihp), \h\^ - + h\. 

To handle the fractional power of the grid operator A, let us consider the eigenvalue 
problem 

We have 

6^A\< 4 < ... < M = (Ni - 1)(N2 - 1), 

where eigenfunctions ip,n, ||</^mll = - 1.2, M. form a basis in H. Therefore 

M 

y = 

m=l 

For the fractional power of the operator A, we have 

M 

A^y = ^Cv, (PmXAiTifm- 

m=l 

Using the above approximations, we arrive from (|4|i at the discrete problem 

A“w = /. (6) 

An approximate solution is sought as a solution of an auxiliary pseudo-time evolu¬ 
tionary problem lE3] . Assume that 

y(t) = (06f(t(A - 661) + 66I)-°‘y(Q), 


with 0 < 0 < 1. Therefore 

y(l) = (0d)“A-“y(O) 

and then w - y(l). The function y{t) satisfies the evolutionary equation 

(tD + 061)—+ aDy ^ 0, 0 < f < 1, (7) 

dt 

where 

D^A- 061. 

By dS, we get 

D^D* >(\-0)6I>O. (8) 
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We supplement (|7]i with the initial condition 


y(0) = (0<5)-“/. 


(9) 


The solution of equation (HJl can be defined as the solution of the Cauchy problem O- 
(|9]l at the final pseudo-time moment t - 1. In 0221] . the case with 0=1 was considered. 

For the solution of the problem Q, (01, it is possible to obtain various a priori 
estimates. Elementary estimates have the form 


IWOlIc < IWO)||g, (10) 

where, e.g., G - I, D. To get (fTOl i for G - D, multiply scalarly equation (0 by dyidt. 
If G = /, then equation 0 is multiplied by ay + tdy/dt. 

To solve numerically the problem 0, 0, we use the simplest implicit two-level 
scheme. Let t be a step of a uniform grid in time such that y" = y(f"), f - nr, 
n - 0,1,..., W Nt - 1. Let us approximate equation 0 by the implicit two-level 
scheme 

«+l _ n 

{r^”'>D + e5r)- - —+aDy°'^"^ n = 0,1,...,A^-1, (11) 

T 

/ = (0d)-“/. (12) 

We use the notation 


rW = -H (1 - 0 -)f", -H (1 - cr)/. 


Lor cr = 0.5, the difference scheme (fTTT i. (fT2l) approximates the problem 0, 0 with 
the second order by t, whereas for other values of cr, we have only the first order. 

We have 


y 


,o-(n) _ 


^ 2 


+ -(y 


1+1 


•/)■ 


Theorem 1. For cr > 0.5 and 0 < 0 < I, the difference scheme Alii . ( 172!) is uncon¬ 
ditionally stable with respect to the initial data. The approximate solution satisfies the 
estimate 


II/"‘|Ig<II/|Ig, « = 0,l,...,iV-l, (13) 


with G - I, D. 


Proof. Rewrite equation (fTTT i as 


H- 061 H- 



2L^,^Z7(y-./) = 0. 

T 2 


Multiplying scalarly this equation by - y", for cr > 0.5, we get 
ll/^‘lb < ll/llfl, n = 0,l,...,iV-l. 

This inequality ensures the estimate (fTsT l for G - D. 
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In a similar way, we consider the case with G - I. Rewrite equation (fTTT l in the 
following form: 


r \ T 


= 0 . 


Multiplying scalarly it by 

in view of ®, we arrive at 


+ fwy_ —L 

T 


< 0 . 


If cr > 0.5, then 
Thus, we obtain (fT3l l for G - I. 


||/"'|| <11/11, n = 0,l,...,iV-l. 


4. Splitting schemes 

The numerical implementation of the difference scheme (fTTT l. (fT^ involves the 
solution of standard elliptic boundary value problems 

+ ao-T)D/+‘ + 0(5y"+' = / 
with the given / for n = 0,1,- 1. 

The inversion of the operator 961 + + Q'crT)(A - 961) may be enough difficult. 

Thus, it seems natural to construct difference schemes for unsteady problems such that 
they will be unconditionally stable, but at the same time, their implementation would 
be considerably simpler. The most interesting results have been obtained taking into 
account a special structure of the problem operator A. 

We dehne a class of the problems (|3l, (|^, where the operator A has the following 
/7-component additive representation: 


p 

(14) 

1=1 

Assume that the operators A,-, i = 1,2,.... /? are simpler than A. We organize compu¬ 
tations in such a way that the transition to a new time level in solving the problem (|7]l, 
(|9]l is not more complicated than the solution of the p problems for the individual oper¬ 
ator terms: A,, i = 1,2,.... /7. Such difference schemes are called additive difference 
schemes 0. 

A decomposition into operator terms in the additive representation (fTTl i may have a 
different nature. In particular, additive difference schemes are used to solve numerically 
multidimensional transient problems of mathematical physics, where one-dimensional 
problems are the most simple ones. On the basis of splitting with respect to the spa¬ 
tial variables, we construct locally one-dimensional schemes. In using computational 
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algorithms of domain decomposition, which focus on modern parallel computers, the 
original problem is divided into several subproblems, and each of them is solved in its 
own subdomain on its individual processor. 

To solve numerically the evolutionary problem (|7]i-(|9]l, it seems reasonable to apply 
the following additive representation of D: 

p 

D = jU/>0, / = 1,2,...,/7, (15) 

1=1 

with pairwise non-permutable operators D/, i - 1 , 2, ..., p: 

DiDj DjDi, i + j, i,j ^ 1,2,.. .,p. 

For the splitting (fT4l l. we put 

p 

Di^Ai-Xil, i^l,2,...,p, 

1=1 

with an appropriate choice of the constants Xi, i - 1,2,..., p. In particular, if 

A, > 6il, 6i > 0, 


then we may take 


Xi-ddi, i-1,2. 


.,p. 


The problem (|7|i, (Ell, (ffSl) does not allow the use of the standard splitting schemes 
because here we have a splitting of the operator at the time derivative. In the work 
1 1811 , we proposed and investigated vector additive-operator schemes with a splitting of 
the operator at the time derivative into the sum of positive dehnite self-adjoint opera¬ 
tors. Here we consider a class of problems with a splitting of both the leading operator 
of the problem and the operator at the time derivative. Such problems are typical in 
studying boundary value problems for pseudo-parabolic equations. Some vector split¬ 
ting schemes for pseudo-parabolic equations with constant operators are considered in 

B . 

For the Cauchy problem (|7|i, (El. (fTSl) . we specify the vector j = {y\,y 2 , ...,yp}. 
Each individual component is determined as the solution of the similar problems 


+ + = o<f<i, 


dy 


dt 


J=i 


dt 


(16) 


J=i 


y,(0) = (05)-“/, i^l,2,...,p. 


(17) 


Here is an elementary coordinate-wise estimate for the solution stability. Subtract¬ 
ing one equation from another, we get 


dyi dyi^i ■ ^ n 

— - i = 2,3,...,p. 

dt dt 
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In view of the initial conditions (fTTl i. we have 




For the individual component y,, we obtain the same equation as for y: 



0 < f < 1, i — 1,2,..., p. 


Thus (see (fTOl i). a priori estimates 


\\yi(mG<m~'^\\f\\c, i^i,2,...,p 


(18) 


are satisfied. Because of this, we have 

yM^yit), 0<t<l, i^\,2,...,p. 

and so, we can treat any component of the vector y{t) as the solution of the original 
problem (|7]i, (|9]l. 

For the vector evolutionary problem (fThl) . (fTTl) . it is easy to obtain (see, e.g., lEl]) 
a priori estimates for the vector y. Let us introduce the Hilbert space H - with the 
scalar product 


p 



Rewrite equation (fThl) as 



0 < f < 1, a - 1,2, ...,p. 


This allows us to write the system of equations in the vector form 



(19) 


Operator matrices B and A seem like this 


B = {Bij], Bij = dijOdDi + tDiDj, 

A = {Aij}, Aij = aDiDj, i,j =1,2, ...,p. 


( 20 ) 


where 6ij is the Kronecker delta. The equation (fT9l l is supplemented with the initial 
condition 


T(0)-/, 


( 21 ) 


where y° = {yi( 0 ),y 2 ( 0 ), ...,yp(0)}. 

The principal advantage of using the formulation (fT^ results from the fact that in 
H, we have 


B = B* > 0, A = A* > 0. 
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Let us consider an a priori estimate for the solution of the vector problem (fT9]l- (l2n i. 
which, on the one hand, is more complicated than (fTSl) . and on the other hand, it will 
serve us as a guide for constructing operator-difference schemes. 

For B, we introduce the representation 


B^C+-A, 

a 


where 

C = [Cij], Cij = 5ij05Di, i, 7 = 1,2,..., p. 
Rewrite equation (fT9l l in the form 

dt \a dt 


By (EOll, we get 


(Ay,y) = a 


and so we have A > 0. Multiply 


\j=i j=i 
scalarly in H by 

t dy 


a dt 


+y- 


This gives 




In view of C > 0, from (l2Tt . we get that the estimate 


0||2 


llTlIc ^ llrll: 

holds. Taking into account the representation for C, we have 


( 22 ) 


(23) 


(24) 


\\l = 06 y (D,y,-,y,) . 


1=1 


Thus, the estimate (l24l i along with ( fTST i may be treated as a vector analogue of the 
estimate ( fTOb . In view of (fTSl l, the estimate (l24l) ensures stability of each individual 
components of the vector y(f). 

To construct splitting schemes for solving the problem (|7]i, (|9]l, we apply common 
schemes with weights and the vector problem ( fT9] l, (|2T]) . In fact, we construct split¬ 
ting schemes for the system of unsteady equations (fT6b , ( fT7b , which are coupled via 
time derivatives. For such problems, the construction of unconditionally stable split¬ 
ting schemes can be caiTied out on the basis of the triangular splitting of the operator 
matrices B and A or the separation of the diagonal part of these operator matrices 
SSIHlIil. In some cases, it is reasonable to combine these approaches, namely, to 
conduct the triangular splitting of the operator matrix B and to separate the diagonal of 
the operator matrix A, or vice versa. 
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Let us construct additive operator-difference schemes using the splitting of operator 
A with separation of the diagonal part. In this case, we obtain 


A = Ao+Au Ao = Amg(Au,A22,...,App). 
In additive representation ( l25l l, we have 



(An 

0 ■ 

0 i 


f ® 

An ■ 

Ayp' 

Ao = 

0 

A 22 ■ 

0 

, Ai = 

A 21 

0 ■ 

■ A2p 


. 0 

0 ■ 

App, 


'.^pl 

Ap2 ■ 

■ 0 . 


In view of 


' ({ P 


= Tj^jyj 


\j=i 


j=i 


/ \V /=1 

P 

< ^ 

/=! 


2 \ 

,1 


i' y 

p Yj ^)"pY i^jyp'^jyj) ’ 


we get 


A < pAo. 

We can consider problem ( fT9b . (|2TI) . (1251) under the additional assumption; 

Ao > -A > 0. 


(25) 


(26) 


To solve the vector problem, we apply the following scheme with weights cri, cr 2 : 




+ — Ao (Ti 


yn+\ _yn 


a 


f v" - 
a T 


+ (1 - CTl) 


yH _yn-\ 


(27) 


-I- A(){a-2y’'^^ + (1 - cr2)y") H- Aij" = 0, 

n = 1,2, 1. 


In the coordinate-wise formulation, the scheme (l27l i corresponds to applying the scheme 


n+l _ n 

QSil -£L + fDi 


( ^ 

y i y i , /1 \ ? y i 

0-1- -H (1 - O-i)- 




yj-yr 


+ aDi(o-2yT^ -H (1 - o-2)y") a ^ DjY] = 0, 

i*j=\ 

n - 0,\, ■■■,N - 

to the Cauchy problem (fThl l. (fTTl l. In contrast to (fTTl l. the scheme (l25l l is three-level 
and it has the weighting parameter cri in the approximation of the time derivative. 
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Computational implementation of dZTl l involves the solution of grid problems 

{0SI+ {a-if + a-2Ta)Di)yi^^ 


for the transition from the time level f" to the new level . 

Using the notation (l20l i. (l25T l. rewrite the operator-difference scheme (l27l i in the 
canonical form of three-level schemes lll5lll6|] : 


«+l _ «-l n + l _ 2yK + 

G"- - - — -H R' - -^- - — -H Ay" = 0. 

2t r 


(28) 


In view of 


y'+l __y" J«+1 J«+1 _ 2^"+' 


2t 


ytl _ ytl—l ~ _1_ 11^-1 


2t 

+y^ 


2t 


2t 


for G and /?, we have 


G" = C +-A + cr2TAo, /?" = Jc + (Tit-Ao - J-A + frz^Ao- (29) 

a 2 a I a 2 


It is essential that the operators G and R are variable, namely, they depend on time. 
When considering the three-level schemes, we have complicated norms, and this un¬ 
steadiness of the operators makes practically impossible to obtain global estimates for 
stability. For this reason, here we formulate a more particular result. 


Theorem 2. If 2cr\ > p and 2cr2 > p, then for the solution of the explicit-implicit 
scheme (O, (O, the following estimate with respect to the initial data 


yH+l _^_yn 

2 

ytl+X _yn 

2 

2 

A 

T 

R 


yn +^n-l 

2 

yn _yn-l 

2 

A 

T 


(30) 


holds. 


Proof. Rewrite (l28T l as 


Introduce 




2t 


n-l 


4 / 

n+l 2y« H- 

-H A- - - - - - = 0. 


1 , y«-y«-l 

v” = -(y"-t-y"-'), W"^- - - - 

2 T 


( 31 ) 
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and rewrite OTI) in the form 


G" 


M'"+‘ + W" 
2 



+ iA(v"+' 


+ v'') = 0. 


Multiplying scalarly (l32l i by 


(32) 


2(v«+' - v") = + tr"), 


we obtain the equality 


+tr'')+ - — a|(h’"+’ +H'" 


(33) 


•(A(v''+' +v"),v"+‘ -v") = 0. 


Taking into account the positive definiteness of the operator G" and the self-adjointness 
of the operators /?", A, from (1^ . we get 

H- (A(v''+', v''+') - (A(v'', v") < 0. 


From this inequality, we obtain dSOl l under the condition 


By (Hg, (iig, we get 



R' - ^A = Ic + ^-(2(riAo - A) + ^-{2(T2Aq - A) 

4 2 2a 4 

> + “(2cri - p)Ao + ^(2o-2 - p)Aq > 0 

2 2a 4 

for 2cri - p >0 and 2 cr 2 - p >0- 


(34) 


5. Numerical experiments 

For simplicity, we restrict ourselves to the case, where ^ is the Laplace operator, 

i.e., 

k(x) = 1, c(jc) = 0 

in equation ([TJ. Under these assumptions, the solution of the spectral problem for 
the differential and discrete Laplace operator is well known and so, it is possible to 
construct the exact solutions. In particular, for the constant 6 in the inequality (|5]l, we 
have 

di =-^sin^:^, k=l,2. (35) 

hi 2Nk 
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Table 1: The eiTor of the two-level scheme with weights for cr = 1 


N 

5 

10 

20 

40 

80 

e(9- 

= 1) 

0.0120292 

0.0066811 

0.0035420 

0.0018307 

0.0009359 

exiO 

= 1) 

0.1362146 

0.0756551 

0.0401086 

0.0207297 

0.0105974 

e(0 = 

0.5) 

0.0255692 

0.0148685 

0.0081424 

0.0042856 

0.0022062 

sa(9 = 

: 0.5) 

0.2869894 

0.1671826 

0.0916372 

0.0482531 

0.0248449 




Table 2: The en'or 

of the two-level scheme with weights for cr = 0.5 


N 

5 

10 

20 

40 

80 

e(0 = 

= 1) 

0.0014351 

0.0004130 

0.0001236 

0.0000484 

0.0000296 

eA(0 ^ 

= 1) 

0.0162505 

0.0046760 

0.0013988 

0.0005462 

0.0003309 

e(0 = 

0.5) 

0.0026070 

0.0008303 

0.0002392 

0.0000720 

0.0000287 

sa (6 = 

: 0.5) 

0.0295158 

0.0094001 

0.0027076 

0.0008137 

0.0003206 


Let the right-hand side is given as 

f{x) - sin(7rxi) sin(7rx2) -H sin(37rxi) sin(27rx2). 

In this case, the exact solution of equation (|4li has the form 

u{x) - sin(7rxi) sin(7rx2) + sin(3;7rxi) sin(27rx2), 
vi = n^, V2 - 

The error of the approximate solution was evaluated in norms H and Ha : 

s^\\y-u\\, eA = ||y-M|U- 

We use the uniform grid in space with - N 2 - 100 and the grid in time with N -20 
as the basic one and the parameter a is equal to 0.5, unless otherwise stated. 

We start with numerical results for the problem (|6]l obtained using (fTTTi . (fT2li . The 
accuracy of the scheme with cr = 1 is presented in Table[T] Similar data for the scheme 
(HU, (IT 2 I with cr = 0.5 are shown in Table |2] It is clear that the scheme of second- 
order accuracy (cr = 0.5) provides good accuracy for a small number of time steps and 
demonstrates a weak dependence on errors in specifying spectrum boundaries for the 
problem operator (compare data at the critical value 0=1 and 6 - 0.5). For the first- 
order scheme (cr = 1), the accuracy is signihcantly lower and influence of 9 is more 
strongly expressed. 

The influence of a spacial grid size on the accuracy of a numerical solution for 
cr = 1.0 and cr = 0.5 is presented in Table[3]for 0 = 1, = 80. In the case of cr = 1.0, 

there is no decreasing of the error with refining spatial grids, i.e., the main part of error 
results from the approximation in time. If cr = 0.5, then error decreases with refining 
spatial grids. Table |4] demonstrates numerical results obtained with various values of 
the parameter a on the grid = A ^2 = 100, N - 80. 
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Table 3: Error for various grids in space 


II 

25 

50 

100 

200 

400 

e(o- = 1) 
sa{o- = 1) 
e(cr = 0.5) 
eA(cr = 0.5) 

0.0010624 

0.0119476 

0.0002321 

0.0025009 

0.0009504 

0.0107503 

0.0000600 

0.0006511 

0.0009359 

0.0105974 

0.0000172 

0.0001888 

0.0009387 

0.0106319 

0.0000066 

0.0000737 

0.0009426 

0.0106768 

0.0000400 

0.0000452 


a 

Table 4: Error for 

0.1 0.3 

various values of a 

0.5 

0.7 

0.9 

e(o- = 1) 
eaCct = 1) 
e(cr = 0.5) 
eA(cr = 0.5) 

0.0009628 

0.0109020 

0.0000925 

0.0002772 

0.0012876 

0.0145807 

0.0000277 

0.0003097 

0.0009359 

0.0105974 

0.0000172 

0.0001888 

0.0005621 

0.0063651 

0.0000090 

0.0000950 

0.0003062 

0.0034674 

0.0000045 

0.0000433 


Possibilities of using splitting schemes for solving the problem (|6ll are shown on 
the splitting with respect to spatial variables, where p = 2 in (fT4l) and one-dimensional 
operators Ai and A 2 are defined as it is given in Section[3 In our model problem with 
the Laplace operator, we have 

Ai>6il, 6i>Q, /=1,2 

with 6i, i - 1,2 defined according to (ILST i. For Di, i - 1,2 (see (flSl l). we put 


Di = Ai -Xil, Xi = i = 1,2. 

In our case of the two-component splitting, the vector scheme (l27l l takes the form 


yT^ - yl 

esil -ii + fDi 


CTl- + (1 - CTl)- 


+ fD2 


yl-yl^ 


+ aDi(o-2y"^' + (1 - cr2)y") -h aD2yl = 0, 


yV^ ~y 2 yl ~y"\ ^ 

OS-il -^ -H fD2 


,,n+l _ \ 

"2 _,_n ^"2 yi 

CTl- H-(l - CTl)- 


-I- aDiy’l -1- aD2{a-2f2*^ + (1 - caK) = 0- 


For three-level schemes, we need to calculate separately an approximate solution at 
the first level using a two-level scheme. Taking into account that we use the scheme 
with the first-order approximation in time, the simplest explicit scheme can be applied. 
Then 

-yO 

QSil —£1 + aDy° = 0, / = 1,2. 

T 

For our two-component splitting with p = 2, we restrict ourselves (see Theorem|2| to 
the case cri = 0-2 = 1. 


15 

















Table 5: The error of the splitting scheme with 6 = \ 


N 

5 

10 

20 

40 

80 

g(i) 

0.0084773 

0.0032158 

0.0012451 

0.0005118 

0.0002210 


0.0959950 

0.0364151 

0.0140991 

0.0057948 

0.0025025 

e(2) 

0.0251082 

0.0080231 

0.0029405 

0.0012031 

0.0005301 


0.2843178 

0.0908516 

0.0332978 

0.0136229 

0.0060020 


N 

Table 6: The error of the splitting scheme with 0 = 0.5 

5 10 20 40 

80 

g(i) 

0.0063711 

0.0045436 

0.0019529 

0.0007352 

0.0002820 

^A 

0.0682378 

0.0503844 

0.0215774 

0.0079984 

0.0029949 

e(2) 

0.1398605 

0.0399402 

0.0130816 

0.0049608 

0.0020928 


1.5835471 

0.4521733 

0.1480800 

0.0561468 

0.0236847 


We can consider the error of individual components of the vector splitting scheme. 
Let 

= Wji - u\l ^ 2. 

In studying the splitting scheme, emphasis is on a dependence of error on a time step. 
We present the results of calculations by the splitting scheme obtained on the basic grid 
in space N\ - N 2 — 100 for a - 0.5. Table|5]shows the accuracy of individual compo¬ 
nents of the vector splitting scheme at the critical value of the parameter 0=1. Similar 
data for 0 - 0.5 are summarized in Table|6] Predictions demonstrate convergence with 
refining grids in time. The accuracy of the splitting schemes is comparable with the 
accuracy of the implicit scheme without splitting (fTTT i. (fTSl i at cr = 1. 
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